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A pseudo-Newtonian Hill problem based on the Paczyriski-Wiita pseudo-Newtonian po- 
tential that reproduces general relativistic effects is presented and compared with the usual 
Newtonian Hill problem. Poincare maps, Lyapunov exponents and fractal escape tech- 
niques are employed to study bounded and unbounded orbits. In particular we consider 
^ ■ the systems composed by Sun, Earth and Moon and composed by the Milky Way, the M2 

Tj- [ cluster and a star. We find that some pseudo-Newtonian systems - including the M2 system 

I - are more stable than their Newtonian equivalent. 
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1 Introduction 
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■ In the 19th century G.W. Hill [1] presented an approximation of the Moon-Earth- 

Sun system, where the movement of the Moon around the Earth was just perturbed 
by a distant Sun. This approximation became known as the Hill Problem. Nowa- 
days, this approximation is still applied in solar system models where bodies in 
nearly circular orbits are perturbed by other far away massive bodies. Also, the Hill 
problem is very useful in the study of the stellar dynamics, e.g., consider a star in 
a star cluster this last orbiting around a galaxy. The star, the cluster, and the galaxy 
can be considered as the Moon, the Earth, and the Sun, respectively. Despite the 
fact that the potentials of the cluster and the galaxy are far from being the poten- 
tials of a point masses and their orbits are far from being circular the Hill problem 
can be taken as a first approximation and can easily accommodate necessary mod- 
ifications, see for instance, Heggie [2]. 
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The Hill Problem was first formulated in the realm of Newtonian dynamics. How- 
ever, there exists extreme cases, like very massive bodies, black holes, and systems 
at great velocities, etc. wherein the Newtonian mechanics is no longer valid and 
relativistic corrections or a fully general relativistic approach is needed. The Hill 
problem was proved to be non-integrable by Meletlidou et al. [3], and is chaotic, as 
shown by Simo and Stuchi [4]. 

The aim of this paper is first to study a pseudo-Newtonian Hill problem obtained 
replacing the Newtonian potential by the Paczyiiski-Wiita [5] potential that repro- 
duces in a rather simple way some aspects of the general relativistic dynamics. 
This is can be considered as a "zeroth order" approach to a general relativistic Hill 
problem. A more rigorous approach should consider a schema of systematic ap- 
proximations within the realm of Einstein theory of gravitation like the schema of 
post-Newtonian approximations. Due to the mathematical complexity of the post- 
Newtonian approximations in this case it is worth to begin with the simpler, but not 
rigorous, treatment of the Hill problem based on the pseudo-Newtonian potential 
above mentioned. The results obtained will be used as a starting point for a more 
complete treatment of the same problem in a future work. 

We shall compare the stability of orbits of the third body in the pseudo-Newtonian 
general relativistic simulation with the equivalent classical Newtonian system. In 
particular we shall study the Sun-Earth-Moon system and the three body problem 
associated with the Milky Way, the M2 cluster and a star. 

In the Section 2 we review the Newtonian Hill problem, with its only integral of mo- 
tion, the Jacobi integral. In Section 3 we present the pseudo-Newtonian Hill prob- 
lem obtained obtained replacing the Newtonian potential by the Paczyhski-Wiita 
potential. In the Section 4 we study the stability of the third body (Moon) orbits. 
The techniques used are: Poincare sections, fractal escapes and fractal dimensions, 
and Lyapunov exponents. Finally, in Section 5, we present our conclusions. 



2 Hill's Equations of Motion 

The Hill's lunar problem is a special case of the circular, planar restricted three- 
body problem, as described by Murray [6] and Arnold [7]. In this problem there is 
a system of two massive bodies with masses mi and m2, m2 << mi , in circular orbits 
around their center of mass and a third massless body moving under influence of 
this system without perturbating it. We can choose units of mass such that G(mi -i- 
mi) = 1. In this way, we can take the masses of the two massive bodies respectively 
as 1 - yu and /i, where ji = Gm2 (note that in this units m2 < mi implies yu < 
1/2) . If we consider the center of mass of the two massive bodies as the origin of 
the coordinate system. The massive bodies with masses // and 1 - yu are situated 
respectively at the positions (1 - ii)R and jiR of a rotating x-axis (see Fig.l), where 
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Fig. 1. The planar, circular restricted three-body problem, on a inertial coordinate system 
if) and on a rotating coordinate system {x, y). Note that on the rotating coordinate system 
the positions of the two massive bodies are fixed. 

R is the distance between the massive bodies. We have 



^ ^ | G(mi +m2) j 



for this distance, where oj is the angular velocity of the rotating frame. By using the 
unit of mass above defined and a unit of time such that co = 1, we have R = 1. In 
these units, the equations of the motion for the third body read 



x = 2y + X - —, (2) 
ox 

y = -lx + y-^, (3) 
dy 

where 

y = -^-^ (4) 

n r2 



with 
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n = ^(x + idf+y^, (5) 
r2 = ^[x-(l-fi)f+y^. (6) 
To obtain Hill's equations, we perform the transformation 

X ^ (I - 1^) + M^'h, (7) 
y^M^'y- (8) 

This transformation corresponds to transfer the origin of the coordinate system to 
the body of mass // and to consider the motion in a disk of radius ju^'^ centred in 
this body. By considering /u small and taken in expansion terms up to OQj.^^^), we 
obtain the Hill's equations 



x = 2y + 3x-^, (9) 
y = -2x-^, (10) 

where r = -yjx^ + y'^. These equations can also be written as 



x-ly = -^, (11) 
ox 

dU 

y + 2x = -—, (12) 
dy 

where U is the Hill potential 
3 . 1 

U = --x'--. (13) 
2 r 



The Hill equations can also be obtained from the Hamiltonian 

h = \[pI + pI) '1 + (yp. - xpy) + \[y^- 2x^) , 



(14) 



where /?^ - x - y and py = j + ;c are the generalized momenta. H is the Jacobi 
Integral, and can be written in the form 

Cj = H = Ux'+f)---\x\ (15) 
2 ^ ' r 2 



Meletlidou et al. [3] proved that Cj = H is the only motion integral of the Hill 
equations, therefore the Hill problem is non-integrable. 
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3 A pseudo-Newtonian Hill problem 



The restricted three-body problem with general relativistic corrections (post-Newtonian 
corrections) is unknown, and as in the simpler case of the two body problem it 
should lead to equations that are not simple for a first approach. By this reason we 
shall simulate general relativistic effects via the pseudo-Newtonian potential intro- 
duced by Paczyhski-Wiita [5] to the study of accretion disks in black holes. Also 
this potential has been used to study black holes with halos by Gueron and Lete- 
lier [8]). Other pseudo-Newtonian models can be found in literature, e.g., the one 
studied by Semerak and Karas [9] to describe rotating black holes. The potential 
introduced by Paczyhski-Wiita [5] changes the usual GM/r Newtonian potential 
by GM/(r - rs), where is the Scharzschild radius, as in the General Relativ- 
ity, rs = IGMjc^. This potential exactly reproduces the marginally bound circular 
orbit, r„^ = Irs, and the last stable circular orbit, r^s = 3rs, and yields efficiency 
factors in close agreement with the Schwarzschild solution. In particular, in the sys- 
tem under consideration, the third body moves in the vicinity of the second body 
with mass jx. The first body with mass 1 - is far away from the second and third 
bodies. Then we change the potential of the second body by yu/(r - rsi), and we 
kept the potential of the first body, (1 - lJ)l{r - rsi) » (1 -fi)/r, where rji and rs2 
are the Scharzschild radius of the first and second body, respectively. The modified 
Hill equations are 



x-2y = 3x-— — ■ = — — , (16) 



r(r-r*)2 dx 
r(r - r*g )^ dy 



y dUp 

y + 2x=-—^-— = -—^, (17) 



where Up is the modified Hill potential 

Up = —x" - (18) 
2 r-r; 



and = //" ' r^a, in units such that the separation between the two massive bodies 
is = 1, like in the Newtonian case. The modified Jacobi constant, in this case, is 
given by 



From the condition that the horizons of the massive bodies do not intersect we have 
an upper limit for r^. This implies r^i -i- rs2 < R- From the condition < 1/2 we 
find (after some algebra) that < 2"^^^ w 0.63. 
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4 Stability of orbits 

We compare the Newtonian system (r* = 0) and pseudo-Newtonian systems with 
parameters varying from 10"^^ to 10"^. The value = 4 x 10"^^ corresponds to 
the Sun-Earth-Moon system, and = 5 x 10"^'' to a system involving the Milky 
Way, the M2 cluster and a typical star. The others represent systems with a very 
massive second body. 

a) Poincare sections 

The orbits in the Newtonian as well in the pseudo-Newtonian Hill problem are the 
solutions of a four dimensional dynamical system with variables {x,y, x,y). Since 
we have an integral of motion, Cj, the motion is reduced to a three dimensional 
system, we can take {x,x,y) as independent variables. We shall study surface of 
section (Poincare sections) evaluating the orbits for different values of the Jacobi 
constant and registering the crossings of the hypersurface y = with y > 0. The 
results for Cj = -2.17 are shown on the Fig. 2. We see that some Kolmogorov- 
Amold-Moser (KAM) tori are destroyed, indicating the transition of the system 
from regular to chaotic behaviour. Note that the picture for = Q and for = 

5 X 10"^ are similar. Another remarkable feature of this maps is the presence, for 
rj = 5 X 10""^ case, of island chains. These structures usually give rise to more 
destroyed tori indicating less stability of the system. 

b) Lyapunov exponents 

We shall study the Lyapunov exponents for the systems above described to better 
analyze the orbits stability. We shall use the Lyapunov characteristic number {X) 
that is defined as the double limit 

X = lim 

f — » oo 

where 6o and 6 are the deviation of two nearby orbits at times and t respectively 
(see AUigood et al. [10]). We get the largest A using the technique suggested by 
Benettin et al. [1 1], in particular we use an algorithm due to Wolf et al. [12]. 

The Lyapunov exponents are not absolute, but dependent on the choice of the time 
scale. We recall that we have fixed the time scale by the requirement that oj = I. 
This defines a time unit that is natural to each particular system. In these time 
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log(5/5o) 



(20) 




Fig. 2. Poincare sections for a) the Newtonian system (r^ = 0) , b) pseudo-Newtonian 
system with r* = 5 x 10""', c) r* ^ 5 x 10"^ and d) r* = 5 x lO"*. The figures a) and c) 
are very similar, and in d) it can be seen the formation of island chains around the stable 
region in the middle. 

units, for two different systems with equal Lyapunov exponents measured in inverse 
seconds and different angular frequencies a>i and 0J2 such that oji < a>2, we will 
have that the system with larger frequency will have small Lyapunov exponent in 
the natural units of time, and vice versa. In other words, in this case, the separation 
of nearby orbits of the third body in each revolution of the first body will be larger 
for the system with smaller a>. So we shall consider that, even though the two 
systems have equal Lyapunov exponents measured in inverse seconds the system 
with smaller angular frequency is more unstable that the one with larger oj. 

For = we get A = 0.185 + 0.003. The Lyapunov exponents for = 4 x 
10-^^, 5 x 10"^°, 5 X 10-^ 5 x 10"^ 5 x lO^^ are shown in the Fig. 3. We have a local 
minimum around = 5 x 10"^*^. 

c) Fractal escape and fractal dimension 

The Poincare sections were obtained for values of the Jacobi constant such that 

the systems are bounded. For values larger than Cjhounded = -2.16 the systems 
are unbounded and the third body escapes. In Fig. 4 we compute contour plot for 
the potential U of the Newtonian system. The pseudo-Newtonian system is quite 
similar and will not be presented her. We see two escape routes, x +00 and 
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Fig. 3. Lyapunov exponents for different values of . Note that there is a minimum around 
= 5 X 10-10. 

X -oo. 

For open systems that have more than one route to escape we can apply the fractal 
escape technique used by Moura and Letelier [13] in the study of the classical 
Henon-Heiles problem. In this method the basins of the escape routes are obtained 
for a set of initial conditions. For chaotic systems, we have the existence of fractal 
basin boundaries (FBB) indicating a great instability of the orbits. In our case we 
chose a subset of the accessible phase space at a fixed Jacobi constant, defined by 
a segment \x\ < a, y = b and < 9 < In, where a and b are constants to be chosen 
appropriately and 6 is the angle that defines the direction of the velocity with respect 
the X-axis. Then the trajectories are integrated numerically, we have three diff'erent 
cases: (1) the body escapes to x ^ +00, (2) it escapes to x ^ -00, and (3) the 
particle does not scape during the integration time. We take the integration time 
long enough to be sure that our results do not depend on this time. Fig. 5(a) shows 
the basins obtained for the Newtonian case. The pseudo-Newtonian basins are very 
similar and will not be shown here. In the figures, grey means initial conditions for 
particles that escape to x ^ +00, black means that the particles escape to x ^ -00 
and white are the initial conditions for particles that do not escape. A zoom of a 
portion of 5(a) is shown in Fig. 5(b), we see self-similarity of the basin boundaries, 
a fractal characteristic. 

To show the diference between systems with diff'erent values of the parameter r* , 
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Fig. 4. Contour plot of the potential U for the Newtonian system. For "energies" above 
-2.16 there are two escape routes, x +co and x —> -oo. 



we calculate the dimension of the basin boundaries obtained for the different sys- 
tems. The dimension used is the box-counting dimension that can be easily ob- 
tained, see for instance, Ott [14] and Grebogi et al. [15]. If we displace a deter- 
mined point of a basin to another on a distance 6, the probability that this new 
initial condition does not belong to the same basin of the old one is, for small e, 
P(e) oc e^"'', where D is the dimension of the set (2 in our case) and d is the box 
counting dimension, also called exterior or fractal dimension when not integer. In 
order to calculate this fractal dimension, for several values of e, we displace the x 
coordinate of all the points from one of the basins (x +oo), and count the number 
of points that does not belong to the same basin. Then we compute the the fraction 
of numbers that does not belong to the same basin, P(e). We plot In P(6) in function 
of In 6. The inclination of the straight line gives us D - d. The obtained values for 
the fractal dimension, d, are shown in the Table 1. Some of the values obtained 
are indistinguishable one from each other in the precision used. Due to his high 
instability, we were not able to calculate the fractal dimension of the r * = 4 x 10""^ 
system. This value is probably near 2, the higher value. 
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Fig. 5. (a) Escape basins for the Newtonian system. Grey (black) means initial conditions 

for particles that escape to x ^ +oo ( x ^ -oo) and white are the initial conditions for 
particles that do not escape during the integration time, (b) A zoom of the portion of the 
precedent figure enclosed by the square. 10 



Table 1 



Fractal dimensions for Newtonian and pseudo-Newtonian systems 



Fractal Dimension 



(Newtonian system) 1 .452 ± 0.003 
4x10-12 1.451 ±0.003 



5 X 10-10 



1.452 ±0.003 



5 X 10-^ 



1.459 ±0.003 



5 X 10-' 



1.466 ±0.003 



5 Conclusions 



We can see directly from the Lyapunov exponents and confirmed by the Poincare 
sections that the pseudo-Newtonian systems (a relativistic simulation) are not al- 
ways more unstable than their equivalent Newtonian systems. 

It could be guessed a priori that the pseudo-Newtonian systems should be more 
unstable due to the fact that the Paczynski-Wiita potential introduces a saddle point 
in the dynamical system. However this happens only for > (9/4)^^^ « 1.31. As 
we mention before, due to the approximation used for the Hill problem, we have 
rj < 0.63. Thus in the cases under consideration the number of saddle points is 
the same for both Newtonian and pseudo-Newtonian systems. Another feature is 
the presence of a local minimum for the Lyapunov exponents, indicating that there 
exists a more stable configuration for pseudo-Newtonian systems that is even more 
stable than the equivalent Newtonian configuration. This minimum corresponds to 
a known physical system, the Milky Way-M2-star system. This is an evidence that 
relativistic effects can made a system more stable. This is related to the fact that the 
addition of extra spherical terms in the Newtonian potential can be used to damp 
the influence of the non-spherical perturbation, making the motion more regular, 
see, for instance, Ivanov et al. [16]. This needs to be confirmed by a fully General 
Relativistic treatment of the Hill problem. We hope to comeback to this subject 
soon. 
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